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Numerical simulations of the cosmic star formation history 
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Abstract. The cosmic star formation history in Cold Dark Matter dominated cosmological scenarios is studied 
by means of hydrodynamical numerical simulations. In particular, we explore a low density model with a A-term 
and two high density models with different power spectra, all of them being spatially flat. Our simulations employ 
a fully nonlinear N-body and Eulerian hydrodynamics algorithm with a model for star formation and supernovae 
feedback that depends on two phenomenological parameters determined in agreement with previous papers. We 
find a nearly constant star formation rate beyond 2 = 1, and we discuss which facts may determine the decrease 
in the SFR from z = 1 to the present epoch. The A-term cosmology with realistic parameters for star formation 
and feedback best reproduces the observed star formation history. 
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1. Introduction 

A detailed understanding of galaxy formation remains one 
of the primary goals of modern cosmology. Whereas grav- 
itation is basically responsible for the large scale struc- 
ture of the universe, many physical phenomena, covering 
a broad dynamical range, are involved on galactic scales. 
A reliable treatment of baryonic physics is essential in or- 
der to obtain a consistent picture of galaxy formation and 
evolution, where star formation and feedback to the inter- 
stellar medium (ISM) play a key role. 

The star formation history of the universe, i.e. the 
global star formation rate (SFR) density as a function 
of redshift, is a crucial test for galaxy formation scenar- 
ios. Many measurements of this quantity have been un- 
dertaken during the past few years, and the number and 
accuracy of available observations are still rapidly increas- 
ing. 

Photometric redshift estimates based on Lyman-limit 
systems (Steidel et al. 1996a, 1996b) stimulated the study 
of the rest frame UV continuum flux around 2000 A as a 
preferred wavelength range for objects in the far universe. 
This emission traces the presence of massive (and there- 
fore young) stars and can be related directly to the actual 
star formation rate. It is relatively easy to detect at high 
z as it is redshifted into the optical band. However, older 
stellar populations and AGN can make a significant con- 
tribution as well, leading to an overestimation of the total 



SFR. On the other hand, dust enshrouding star formation 
regions absorb very efficiently the UV light, re-emitting 
it in the far IR. Unfortunately, there is still some uncer- 
tainty to account for dust extinction and its evolution with 
redshift (Nagamine et al. 2001| ). 
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Recombination lines from H II regions can be in prin- 
ciple a more reliable estimator of the instantaneous SFR, 
since the ionizing radiation (A < 912A) comes from more 
massive (younger) stars than the softer UV continuum. 
These optical lines are less dramatically affected by dust 
extinction and can be mapped with higher resolution in 
the local universe, but they require infrared spectroscopy 
at moderate redshifts, which has not been possible until 
very recently due to the faintness of the sources involved 
(Pettini et al. 1995). For H Q , this happens at z ~ 0.5, and 
for forbidden lines, such as [OII]3727A, at z ~ 1.5. There 
is also a considerable degree of uncertainty in the relation 
between recombination lines and star formation activity. 

Finally, the far infrared (FIR) spectrum of a galaxy 
can be roughly separated into two components, namely 
the thermal emission of dust (heated by star formation 
bursts) around A ~ 60 /xm and an infrared cirrus pow- 
ered by the global radiation field, which dominates for 
A > 100 fim. Neglecting cirrus and AGN contamination, 
the FIR luminosity would also be an excellent tracer, but 
instrumentation is not as developed in the IR and radio 
bands as it is in the optical, and the conversion factor be- 
tween FIR luminosity and SFR is for the moment rather 
model dependent. 
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Although there is a large scatter among different 
indicators at any given redshift, basically all studies find 
a significant increase (by about one decade) in the cosmic 
SFR density from the present day to z = 1. However, it 
is still a controversial issue (Hopkins ct al. 2001) whether 
it reaches a broad maximum there (Madau et al. 1996; 



is included with an analytical prescription for the trans- 
formation between cosmological models. Sect. || describes 
the code and the different numerical experiments on which 
the present paper is based. Our results are presented in 
Sect. |], where we compare with observations and evaluate 
the separate contribution of cosmology, feedback effects 
and environment. Our main conclusions are summarized 
rcdshifts (]Sawicki et al. 1997| ; [Pascarelle et al. 1998| ; in Sect. [|. 



Gispert et al. 2000) or it declines gradually towards higher 



Steidel et al. 199S) 



From a theoretical point of view, star formation is a 
highly non- linear process, which precludes any simple an- 
alytical treatment. Most efforts towards modeling the star 
formation history of the universe resort to semi-analytical 
methods or numerical simulations to tackle the feedback 
mechanisms that self-regulate the SFR within the hierar- 
chical clustering scenario of struct ure formation. 

Semi-analyti cal models ( Kauffmann et al. 1993 ; 
Cole et al. 1994 ) construct a large sample of Monte 
Carlo realizations of halo merging histories, using 
the Press-Schechter formalism. Hydrodynamics, star 
formation and feedback are implemented through sim- 
ple recipes, whose parameters are fixed in order to 
match the observed properties of real galaxies, specially 
the luminosity function and the Tully-Fisher relation 
( (Bomcrvillc fc Primack 1999| ) . Their high computational 
efficiency allows a fast exploration of the parameter 
space in search of a viable model, but the number of free 
parameters asks for supplementary modeling based on 
more fundamental physics. 

Numerical simulations are intended to solve directly 
the set of differential equations of gravitation and gas 
hydrodynamics, and hence require much fewer model as- 
sumptions than the semi-analytic approach. However, de- 
spite the staggering advances in computer technology and 
numerical algorithms, resolution in large-scale hydrody- 
namical simulations is severely restricted by computa- 
tional and data storage requirements. Although gas can be 
treated self-consistently, some phenomenological recipes 
are still required in order to model star formation and 
feedback. 

The cosmic SFR density is a direct outcome of semi- 
analytical methods ( [gomerville et al. 2001 ) and hydrody- 
namical simulations ( Nagamine et al. 2001 ). Both of them 
predict a significant star formation activity at high z, fa- 
voring the hypothesis of a gradual evolution of the comov- 
ing SFR as a function of time. 

The present analysis is based on a series of nu- 
merical simulations, within different cosmological mod- 
els. In previous papers, we have studied basic observ- 
able properties of simulated galaxy samples, such as the 
luminosity function and color dependence on environ- 
ment (Elizondo et al. 1999a), and the Tully-Fisher rela- 
tion ( [Elizondo et al. 1999h| ) In this paper, we will focus 



on the evolution of the SFR on a cosmological scale, com- 
paring our results with recent observations. 

The following section is devoted to the observational 
estimation of the cosmic SFR. The standard procedure 
is briefly summarized and a compilation of recent results 



2. Observations 

The first step to obtain an estimate of the SFR density 
from observations is the selection of a complete sample 
of galaxies at a given redshift and wavelength. Some au- 
thors (Steidel et al. 1999) claim that the HDF covers a 
very small area (~ 5 arcmin 2 ), and hence it is not statis- 
tically representative of the whole universe. This cosmic 
variance is difficult to quantify, but the strong clustering 
of high redshift galaxies (and the fact that most estimates 
beyond z = 2 are based on the HDF) seems to indicate 
that it could be an important effect. 

Another crucial issue is the correction for incomplete- 
ness in the construction of the luminosity function, which 
can be done following the V max formalism ( Schmidt 1968 ) 
or resorting to more elaborate Monte Carlo algorithms 
(e.g. Steidel et al. 1999). The comoving luminosity den- 



sity is easily obtained by integration over all magnitudes 
(extrapolating the faint end as a Schechter function), and 
then it can be converted to a SFR density with the help 
of a population synthesis model. 

A compilation of recent results is given in Table |. The 
first column indicates the bibliographic reference, as well 
as the survey from which the galaxy sample was extracted. 
Next columns correspond to the chosen wavelength and 
redshifts of each analysis, and finally the comoving lumi- 
nosity and SFR densities. 

Conversion factors between the last two quantities dif- 
fer significantly from author to author, and in some cases 
only the luminosity density was provided. If there is a SFR 
estimate in the original paper, then it is quoted in Table [l] 
as corrected by dust extinction] else, we follow Kennicutt's 
(1998) prescription for an exponential burst and a Salpeter 
IMF between 0.1 and 100 M Q : 



P* [ M o yr" 



1.4 x 10- 28 iuv [erg s^Hz- 1 ] 
1.4 x 10- 41 L 37 2 7 [erg s^Hz- 1 ] 



(1) 



For dust extinction, we have applied a correction of 
A(1500-2000A)=1.2 mag and A(2880A,O II)=0.625 mag, 
which correspond to factors of 3.02 and 1.78 respectively. 
All modifications to the original data are shown as a foot- 
note on the bibliographic reference in the Table |l|. 

Most observational values given in Table | were cal- 
culated assuming a standard cold dark matter (SCDM) 
cosmology with dimensionless Hubble parameter h = 0.5 
and mean matter density f2 m = 1. Since the conversion 
factor from luminosity to SFR and the correction for dust 
extinction are completely independent of the cosmological 
scenario, the computation of the comoving SFR density in 
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HDF (Pascarelle et al. 1998) 



HDF (Madau et al. 1998 



(Steidel et al. 199S) 



( fTreyer et al. 1998|) d 
CFRS ( [Lilly et al. 199$ d 



HDF (Cowie et al. 199S) 



HDF (Connolly et al. 1997) 



HDF (Sawicki et al. 1997) 



(GalleRO et al. 1995) 



CFRS ([Tresse fc Maddox 1998| ) 
CFRS 



(IClazcbrook et al. 1999, ) 
(|Yan et al. 1999|) 



CFRS (Hammer et al. 1997) 



CFRS ( [Flores et al. 1999j ) 



HDF (IHuRhes et al. 1998|) 



(Haarsma et al. 2000) 
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-1.23 
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-1.51 


-1.59 


1500 A 


2.75 ± 0.75 


26.90 ± 0.15 


-0.95 


-1.02 




4.00 ± 0.50 


26.50 ± 0.20 


-1.35 


-1.43 


1700 A 


3.04 ± 0.35 


27.04 ± 0.07 


-0.81 


-0.89 




4.13 ± 0.35 


26.95 ± 0.10 


-0.90 


-0.98 


2000 A 


0.15 ± 0.15 


26.45 ± 0.15 


-1.40 


-1.35 


2800 A 


0.35 ± 0.15 


26.14 ± 0.07 


-1.71 


-1.68 




0.625 ± 0.125 


26.46 ± 0.08 


-1.39 


-1.39 




0.875 ± 0.125 


26.78 ± 0.15 


-1.07 


-1.10 


2800 A 


0.35 ± 0.15 


25.91 ± 0.12 


-1.69 


-1.66 




0.625 ± 0.125 


26.15 ± 0.12 


-1.45 


-1.45 




0.875 ± 0.125 


26.06 ± 0.12 


-1.54 


-1.57 




1.25 ± 0.25 


26.35 ± 0.10 


-1.25 


-1.30 


2800 A 


0.75 ± 0.25 


26.77 ± 0.15 


-1.08 


-1.10 




1.25 ± 0.25 


26.94 ± 0.15 


-0.91 


-0.96 




1.75 ± 0.25 


26.84 ± 0.15 


-1.26 


-1.32 


3000 A 


0.35 ±0.15 


26.51 ± 0.47 


-1.34 


-1.31 




0.75 ±0.25 


26.74 ± 0.06 


-1.11 


-1.13 




1.50 ± 0.50 


26.93 ± 0.05 


-0.92 


-0.97 




2.50 ± 0.50 


27.28 ± 0.06 


-0.57 


-0.64 




3.50 ± 0.50 


26.91 ± 0.10 


-0.94 


-1.02 


H Q 


0.025 ± 0.025 


39.09 ± 0.04 


-2.01 


-1.89 


H„ 


0.2 ± 0.1 


39.44 ± 0.04 


-1.66 


-1.60 
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0.875 ± 0.125 


40.01 ± 0.15 


-0.91 


-0.94 


H Q 


1.3 ± 0.5 


40.21 ± 0.13 


-0.89 


-0.94 


II 3727 A 


0.35 ± 0.15 


38.63±°;°| 


-1.97 


-1.94 




0.625 ± 0.125 


39.16to.iJ 


-1.44 


-1.44 
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-1.04 


-1.0/ 


15 fim 


0.35 ±0.15 


41.97 ±0.25 


-1.46 


-1.43 




0.625 ±0.125 


42.20 ±0.25 


-1.15 


-1.15 




0.875 ±0.125 


42.53 ±0.25 


-0.82 


-0.85 


60 fim 


3.0 ± 1.0 


< 42.26 


-0.98 


-1.05 


1.4 GHz 


0.28 ±0.12 


26.75+c^ 


-1.17 


-1.13 




0.46 ±0.05 


27.03l ; 2 2 t 


-0.89 


-0.87 




0.60 ±0.05 


97 1 O+0-16 


-0.80 


-0.80 




0.81 ±0.08 


27.39«; 1 1 3 8 


-0.53 


-0.55 




1.60 ±0.64 


27.54«; 1 1 3 8 


-0.38 


-0.44 



Table 1. Observational estimates of the cosmic SFR density at different epochs. /9l refers to the luminosity density at 
the appropriate wavelength, expressed in erg s _1 Hz _1 Mpc -3 . represents the comoving SFR density (in M yr^ 1 
Mpc^ 3 ) for our SCDM and ACDM cosmologies. d Original data has been corrected for dust extinction. 



a ACDM (h = 0.7) model only involves the transformation 
of the luminosity densities. 

In a flat universe, the volume enclosed by a solid angle 
Auj between z — Az and z + Az is 



4tt 



V(z, Az) = —Auj (d 3 m (z + Az) - <(z - Az)) 



(2) 



dh{z) = (l±z)d m (z). Then, assuming that all galaxies are 
located at the centers of the interval in z, the luminosity 
density (and hence the total SFR density) will be propor- 
tional to 



where d m (z) is the comoving distance in each cosmological 
model, while luminosity scales with luminosity distance as 



P*{z) oc 



L(z) 



dl(z) 



V(z,Az) d^(z + Az) - df n (z - Az) 



(3) 
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Fig. 1. Conversion factor between star formation densities 
in our cosmological models SCDM (h = 0.5, f2 in = 1, 
A = 0) and ACDM (h = 0.7, Q m = 0.35, Q A = 0.65) for 
several redshift intervals Az. Dashed line corresponds to 
the approximation (||). 

To first order in Az, the conversion between SFR den- 
sities in ACDM and variants of the SCDM (h = 0.5) cos- 
mologies can be approximated as 

^ACDM ^ACDM 
pSCDM ~ ^SCDM 



y/Sl m + Sl A (l + z)- 



(4) 



Fig. [l] shows this factor, as well as the exact expres- 
sion for different values of Az, which requires a numerical 
evaluation of the luminosity distance in the ACDM model. 
Although the results are very similar, we have chosen to 
use the exact expression (||) to compute the SFRs in Table 
|l| in order to minimize errors. 

At low redshifts, the data are consistent for all the in- 
dependent tracers of the star formation activity once the 
effect of dust obscuration is taken into account. The steep 
increase in the cosmic SFR density from the present day 
to z ~ 1 is firmly established. However, there is not yet 
a full agreement about whether there is a characteristic 
epoch of star formation in the universe around z <~ 1.5 
or on the contrary that the SFR declines very slowly to- 
wards higher z, as more recent observations seem to point 
out. The results of ongoing deep surveys will be extremely 
helpful in order to reach a definitive conclusion from the 
observational point of view. 

3. Simulations 

3.1. Cosmological Models 

Up to now, the ACDM model has proven to be very suc- 
cessful in describing most of the observational data at 



both low and high redshifts. Moreover, recent observa- 
tions indicate the existence of a positive cosmological con- 
stant QPcrlmutter et al. 199S| ; fVlclchiorri et al. 2000| ). For 
our study we have chosen a spatially flat universe with a 
cosmological constant A/3Hq = Oa = 0.65 and a present- 
day Hubble constant of Ho = lOO/i = 70 km/s/Mpc. In or- 
der to study the influence of cosmology on the global SFR, 
we have investigated three other scenarios: ACDM with 
lower Ha = 0.6, the standard CDM model (dominated by 
dark matter density), and the BSI model with a Broken 
Scale Invariant perturbation spectrum (as predicted by 
double inflation, Gottlobcr et al. 1991 ). The latter model 



is also dark matter-dominated, but it describes a primor- 
dial power spectrum in which fluctuations on small scales 
are significantly suppressed. It is almost the same spec- 
trum as in the rCDM model, in which the primordial per- 
turbation spectrum has been changed due to the hypothet- 
ical decay of massive r neutrinos ( Efstathiou et al. 1992 ) . 
The model parameters are summarized in Table |2[ All 
spectra are COBE-normalized. The normalization is given 
by the linear dark matter mass variance as in a 8/i -1 Mpc 
sphere extrapolated to present time. 

3.2. Numerical code 



The code used to performed the numerical experiments 
consists of a combination of a Particle-Mesh N-body grav- 
ity solver and an Eulerian Piecewise Parabolic Method 
to account for gasdynamics. A detailed modeling of the 
most relevant non-adiabatic baryonic processes (i.e. mul- 
tiphase ISM, cooling, photoionization, star formation and 
supernovae feedbacks) acting at a sub-grid level is also 
i ncorporated. The code is described in detail elsewhere 
( Yepes et al. 1997] ). Here we briefly review its main fea- 
tures regarding star formation and feedback. 

The ISM is modeled as a two phase medium: hot gas 
(T h > 10 4 K) and a cold phase (T c = 10 4 K), which pre- 
tends to mimic the molecular clouds from which stars are 
formed. 

Cold clouds can be formed from the hot gas by ra- 
diative and Compton cooling or by thermal instabil- 
ity. Radiative cooling rates, A r (T, Z), are taken from 
[Sutherland fc Dopita (1993) for two metallicity cases: ei- 
ther primordial composition or solar, depending on the 
star formation history of each computational element. 

Stars are created, from the cold phase, as collisionless 
particles at each timestep. Their mass is proportional to 
the density of cold gas in each cell: 



p* 



(5) 



where f3 is the fraction of massive stars (M > 10M Q ) that 
explode as supernovae. For a Salpeter IMF, (3 = 0.12. For 
the characteristic time of star formation t* we assume a 
constant value of 10 8 yrs. 

Several feedback mechanisms are included: photoion- 
isation by UV light, gas heating and evaporation by su- 
pernova explosions. Massive stars are assumed to create a 
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Experiment 


n r 




fib 


n A 


h 


08 


N 


£box 


Lccll 


Mom 


A 


V 


ACDM1 


1 


0.325 


0.025 


0.65 


0.7 


0.9 


350 3 


30 


85.7 


28.1 


200 


100 


ACDM2 


1 


0.325 


0.025 


0.65 


0.7 


0.9 


270 3 


21.4 


79.3 


22.2 


200 


100 


ACDM3 


1 


0.325 


0.025 


0.65 


0.7 


0.9 


300 3 


12 


40 


2.85 


200 


100 


ACDM4 


9 


0.37 


0.03 


0.6 


0.6 


0.9 


128 3 


5 


39.1 


2.22 


200 


100,500,1000 


ACDM5 


18 


0.324 


0.026 


0.65 


0.7 


0.9 


128 3 


5 


39.1 


2.65 


0,50,200 


100 


SCDM 


18 


0.95 


0.05 





0.5 


1.2 


128 3 


5 


39.1 


3.96 


0,50,200 


100 


BSI 


18 


0.95 


0.05 





0.5 


0.6 


128 3 


5 


39.1 


3.96 


0,50,200 


100 



Table 2. Collection of cosmological models, simulation specifics and feedback parameters. The box length Lbox is 
given in h" 1 Mpc, the cell length L cc \\ is given in h~ x kpc and the mass of dark matter particles is given in 10 6 M Q . 



homogeneously distributed background of ionizing radia- 
tion, and gas should be dense enough to screen it in order 
to be able to cool into the cold phase and form stars (e.g. 



Giroux fc Shapiro 1986 Miicket fc Kates 1997). A simple 



prescription, without resorting to a detailed computation 
of photoionization heating is to impose an overdensity 
threshold T> for star formation to take place. On the other 
hand, supernova explosions can evaporate the cold gas 
phase and increase the temperature of the hot compo- 
nent. In our model, the supernova feedback parameter A 
controls the amount of energy (relative to the supernovae 
mass) that goes into cold cloud evaporation: 



SN 



(Pc) SN = —^A 



0£c 



SN 



(e S N + Ae c ) 



(6) 



(7) 



where the subscripts h and c denote hot and cold gas, 
esN — 4.45 x 10 49 erg Mg 1 is the supernovae specific energy 
(Salpeter IMF) and e c is the specific energy of the cold gas 
phase (T c = 10 4 K). 

3.3. Numerical experiments 

Table || summarizes our set of experiments. In order to 
improve statistics, several realizations were run when the 
number of particles was low enough (N = 128 3 ), changing 
the random seeds used to generate the initial conditions. 
The number of realizations is given by n T in the second 
column. The cosmological model and the normalization u% 
are shown in the next 5 columns, as well as the number of 
particles, box size, cell size (which determines our spatial 
resolution) and mass of each dark matter particle (mass 
resolution). In the last two columns we present the values 
of the supernova feedback parameter A and the overden- 
sity threshold T> required for star formation. 

The comparison of experiments ACDM5, SCDM and 
BSI allows us to determine the effect of cosmology and 
feedback from supernovae on the global SFR history. 
Effects of photoionization are discussed from the results 
of experiment ACDM4, and the dependence on resolution 
and volume can be studied from the whole set of ACDM 
experiments. In ACDM 1 and 2, the volume is high enough 



to include a cluster of galaxies, so we have also investigated 
the behavior of the SFR in different environments. 



3.4. Resolution effects 

In pure A^-body simulations the lack of mass and force 
resolutio n leads to the well known overmergin g problem 
(see e. g. Klypin et al. 1999 ; van Kampcn 2001 ), i. e. sub- 
structures of large halos or halos inside of halos are erased. 

This effect becomes important in galaxy clusters. To 
avoid overmerging high force and mass resolution are nec- 
essary which may be reached later in adaptive mesh and 
multimass N-body/gasdynamical codes. 

Our code is based on a fixed Eulerian mesh, so resolu- 
tion is limited to the cell size. High spatial resolution can 
only be achieved by increasing the number of cells uni- 
formly, or reducing the computational volume. As a trade 
off between the two, we have decided to restrict us on 
medium spatial resolution, in favour of large simulation 
volumes. Resolution effects are discussed in detail in the 
next section. 

The numerical code has been parallelized using 
OpenMP compiler directives. Therefore, it runs very ef- 
ficiently in parallel computers with either Shared or Non 
Uniform Memory Access (NUMA) architectures. The sim- 
ulations have been performed on a variety of machines. 



4. Results 

4.1. Cosmology dependence 

The global star formation histories from simulations of 
the three cosmological scenarios considered are compared 
in Fig. p[ The initial power spectrum normalization de- 
termines the evolution of structures, which represent a 
crucial factor in determining the rate of conversion of gas 
to stars. 

Cosmological scenarios with CDM alone cannot ex- 
plain structure formation on both small and very large 
scales. However, we have included SCDM as a reference 
model in our comparison. Due to its large power at small 
scales, SCDM produces many structures at high redshift. 
Therefore, we get a flat star formation rate even beyond 
z ~ 6, with a plateau about one decade over the observa- 
tions (Fig. ||, middle panel). It does not seem very likely 
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Fig. 2. Evolution of the comoving SFR density in different cosmologies (experiments ACDM5, SCDM and BSI). The 
average value over different realizations is plotted as a thick solid line, while the shaded area represents the la deviation 
from the mean. Dots correspond to the observational data points listed in Table pi 



that this difference can be explained by dust absorption. 
In this model, too many baryons are transformed into stars 
too early with respect to observations, which poses an ad- 
ditional problem for the SCDM model. 

On the contrary, the BSI model leads to a clear max- 
imum in the SFR at z ~ 2, but the absolute value is too 
small and the decrease at higher redshifts is clearly too 
steep (Fig. ||, right). This model has not enough power at 
small scales for producing enough stars at high redshift. 
As a side effect, the cosmic star formation rate density re- 
mains too high at low redshifts, since a large part of cold 
gas remains in high density regions triggering further star 
formation activity. 

A cosmological model with non-zero cosmological con- 
stant offers the best agreement with observations (Fig. g, 
left) . Both the overall trend and the height of the plateau 
fit reasonably well to the observational data points. At 
present, the AC DM model fits bes t all known observa- 
tional data (e.g. Bahcall et al. 1999). Our results confirm 
that this model can also reproduce the observed star for- 
mation history of the universe when dust extinction is 
taken into account, in which case a significant amount of 
star formation is expected to take place at high redshift. 

A further feature (common to all cosmologies) is the 
smooth redshift dependence of the mean SFR and the 
small range of scatter. This is in sharp contrast to the 
evolution of individual baryonic clumps, or the behavior 
of different regions in a single time-step of the simulation. 
There we observe huge bursts of star formation with an 
increase in activity by a factor of 10 and higher, as well as 
a quiescent regime with a star formation rate lower than 
the average. 

In these simulations, the first stars are created around 
z ~ 10, which seems to be not too unrealistic, but it may 
depend on resolution. This was not a primary issue of the 
present study. 

4.2. Feedback dependence 

Fig. H shows the results of experiments ACDM5 and 
ACDM4, from where the effects of supernova feedback and 



photoionization were studied. The same trends described 
in this section were also seen in SCDM and BSI simula- 
tions. 

Supernovae (left panel of Fig. ||) play a key role in the 
slope of the SFR density at low z. Gas heating and evap- 
oration act as a self-regulating mechanism that inhibits 
subsequent gas infall when many new stars are formed. 
This is specially important in the most massive objects, 
which show a much higher star formation activity when 
feedback is not implemented. Supernova explosions can 
also blow away the gas reservoirs of dwarf galaxies. The 
effect on the global star formation rate is more notice- 
able at high redshift, when there are still very few galaxies 
massive enough to retain the heated gas within their grav- 
itational potential wells. However, photoionization (right 
panel) is the dominant effect at very early epochs, prevent- 
ing the gas from cooling and forming stars in low density 
regions (Efstathiou 1992). As we have discussed above, 
our model takes into account this effect through the over- 
density threshold parameter T> for star formation. Only 
in regions denser than T>, molecular clouds will be able to 



-0.5 




No Feedback 
A=200 



6 
z 



Fig. 3. Effects of stellar feedback on the cosmic SFR den- 
sity, a) ACDM5 simulations with supernova feedback pa- 
rameter A — 200 (solid line) within la error (shaded area) 
and A = (dashed line), b) Same as (a) but for ACDM4 
simulations run with different values of the overdensity 
threshold D imposed by photoionization. 
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form, and the onset of cosmic star formation activity is 
thus retarded until more and more objects become dense 
enough to screen the photoionizing background. 

We found a SFR density consistent with observational 
data for T> ~ 100. This is in roughly agreement with re- 
sults from a more detailed theoretical estimate of the same 



parameter ( Miicket fe Kates 1997 ), although in that case 
the actual value of T>, as well as the temperature range 
in which thermal instability is an efficient process, depend 
on the local gas density and the intensity of the ionizing 
flux. 



4.3. Environmental dependence 

In the hierarchical scenario of structure formation, the first 
objects collapsed around the highest overdensity peaks 
and merged to eventually form the giant ellipticals in the 
centers of clusters and the bulges of the most massive iso- 
lated galaxies. These objects contain the vast majority of 
stars older than z = 3, and hence their star formation 
history must be completely different from that of galaxies 
living in less dense environments. 

The experiments ACDM 1 and 2 sample a volume con- 
taining clusters, filaments, and groups of galaxies as well 
as isolated galaxies. Gottlober et al. (2001) have studied 
the merging history of dark matter halos as function of 
halo environment. They found that halos located inside 
clusters have formed earlier than isolated halos of the 
same mass. Moreover, they showed that at higher red- 
shifts (z ~ 1 — 4), progenitors of cluster and group galaxies 
have 3-5 times higher merger rates than isolated galaxies. 
Mergers of galaxies are thought to play a crucial role in the 
evolution of galaxies. In particular, the inflow of material 
may serve as a source of fresh gas and therefore increase 
the star formation rate. 

In Fig. H we plot the contribution to the cosmic SFR 
density of progenitors of objects more massive than 10 13 
Mq at present (clusters, large groups) and less massive 
(isolated galaxies or small groups). The most massive ob- 
jects have their peaks of star formation activity around 
z ~ 3, showing passive evolution since z ~ 2. This is in 
excellent agreement with the merging history found by 
Gottlober et al. (2001)| . The progenitors of these objects 



were responsible for the bulk of star formation at high red- 
shift, and can be very likely associated with the Lyman 
Break Galaxy population. 

Isolated galaxies and galaxies in small groups are 
younger. They are responsible for most of the current stel- 
lar production. In these areas, halos massive enough to 
retain their gas content have a nearly constant SFR (with 
sporadic burst episodes). 

The simulated volume in our numerical experiments is 
too small to be a representative cosmological volume. In 
fact, the larger volume ACDM1 contains only one galaxy 
cluster of approximately 3 x 10 14 M Q whereas the sam- 
ller volume ACDM2 contains four massive groups of about 
5 x 10 13 Mq. Comparing the solid lines in Fig. [|we find 




Fig. 4. Star formation rate densities in models ACDM1 
and ACDM2 due to halos more (solid lines) and less 
(dashed lines) massive than 10 13 M Q at z = 0. 



qualitatively the expected scenario: the peak of star forma- 
tion in the progenitors of the more massive cluster hap- 
pened earlier (about z = 4.5) than for the less massive 
groups (z = 3). 

In the absence of stellar feedback, most of the bary- 
onic matter would have cooled into small mini-galaxies 
( |Cole lgglt |White fc Frenk 199lfc |Blanchard et al. 1992j ). 
Photoionization and supernova explosions inhibit cooling 
and star formation in very low-mass objects, so early star 
formation was biased towards the high overdensity regions 
that will collapse later to form the cores of the first clus- 
ters. 



As the collapse proceeds, gas is heated to tempera- 
tures that make the cooling times too long and the SFR 
drastically drops once the cold gas has been exhausted. 

At the same time, new galaxies will form and pro- 
duce stars. Some of them will remain in rather iso- 
lated areas until present (dashed lines in Fig. [|) and 
some others will fall into the potential wells of clus- 
ters. It seems very likely that tidal interactions trigger 
intense star formation bursts in these infalling galaxies 



(Lavery & Henry 1988), but ram pressure of the ICM 



may also deprive them from their c old gas reservoirs 
( |Gunn fc Gott 1972[ |Quilis et al. 2000|) . This later mech- 
anism is supported by the observational evidence of 
the displacement of the neutral hydrogen disk of post- 
starburst galaxies with respect to the stellar population 
( IValluri fc Jog 1991|) . 

Cluster galaxies would then stop/decay their star for- 
mation activity and become redder, giving rise to the 
[Butcher fc Oemler (1978) effect. At the same time, they 
will experience a morphological transformation, due to 
mergers wit h their neighboring c luster galaxies, from spi- 
rals into SO ( Dressier et al. 1997 ). At present, most o f the 
observed clusters are at z < 0.5 ( Fasano et al. 200C , and 
references therein), so deeper surveys are clearly needed 
in order to check this scenario. 
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Fig. 5. Comoving SFR in our 5 ACDM experiments, all 
with A = 200 and T> — 100. A significant star formation 
activity is found at high z regardless of resolution and 
volume. This is consistent with the full set of observational 
data points (Table |l|) . 



4.4. Resolution and volume dependence 

Star formation and feedback processes are incorporated 
in gasdynamical simulations (either SPH or Eulerian) by 
means of rather simple recipes which pretend to extrapo- 
late the effects of local physics acting on scales of a few pc 
and masses of ~ 10 5 M Q below the resolution limit (often 
at kpc and higher). In our case, part of this extrapolation 
is hidden in the two parameters: The supernova feedback 
parameter A and the overdensity threshold for star for- 
mation T>. Whether this parameterization is a reasonable 
extrapolation of the underlying subgrid physics depends 
on the stability of results against changes in spatial reso- 
lution of the simulations. 

On the other hand, when one tries to estimate cos- 
mological averaged quantities, like the SFR density, it is 
necessary to check for deviations of these quantities due 
to small number statistics. In the case that matters here, 
a volume averaged quantity, the cosmic variance can have 
a relatively large contribution in the particular determi- 
nation of the comoving SFR density. Thus, we have sim- 
ulated different computational volumes with similar reso- 
lutions, depending on the computational facilities. 

To check for these two effects we compare in Fig. |^ our 
results for ACDM simulations with different resolutions 
and volumes. 

The most noticeable feature that we derive from this 
Figure is the steeper slope of the local SFR for larger vol- 
ume simulations. As we increase the volume, larger struc- 
tures are being simulated (i.e. clusters) and as explained 



above, the reduced star formation activity in the cores of 
clusters in recent epochs can explain this behavior. In sim- 
ulations with smaller volumes, most of the galaxies that 
form are either isolated or in small groups, which have a 
completely different star formation activity. Thus, we can 
conclude from our simulations that the effects of cosmic 
variance are quite important in determining the proper 
behavior of the cosmic SFR density from z > 1 to the 
present. 

We do not find significant differences in our estimates 
of the SFR density from simulations with different vol- 
umes and spatial resolutions, but there is a general trend 
to have more SFR density at early times for simulations 
with higher resolution. This is a consequence of the over- 
cooling problem (see Balogh et al. 2001 for a recent re- 
view), since the number of resolved low mass halos in- 
creases in these simulations. A detailed analysis of this 
question would require a more realistic treatment of pho- 
toionization. 

Keeping all these considerations in mind, we should 
be extremely cautious when making quantitative estimates 
of the cosmic SFR density at any given z. However, this 
quantity shows a similar behavior in all ACDM experi- 
ments. This is a very promising result indicating that its 
qualitative evolution (particularly at high redshift) is a ro- 
bust prediction of our simulations. 

No characteristic epoch of cosmic star formation is 
found for a ACDM universe. The comoving SFR density 
increases with lookback time until it reaches a nearly con- 
stant mean plateau beyond z — 2. Stars probably formed 
gradually at high redshift, according to most recent ob- 
servations and taking dust extinction into account. At 
some point around z ~ 1, the global star formation ac- 
tivity would have declined sharply until the present day, 
although individual SFRs of active galaxies are kept ap- 
proximately constant. 

The star density p*(z) shows a similar behavior for 
both CDM and ACDM models. The most remarkable dif- 
ference is found in the present value of p*(0) (see Table |J), 
which can be observationally deduced from the local lu- 
minosity density, assuming a constant mass-to-light ratio 
([Madau et al. 1996|; |Fukugita et al. 1998|): 



P° h s (0) ~7xlO 8 /iM Mpc" 



(8) 



Base d on the fracti on of stars in bulges, ellipticals 
and SOs, Renzini (1998) estimated that about 30% of the 
stellar content of the universe should have already been 
formed by z = 3, regardless of the cosmological model. 
The amount of stars produced in experiments ACDM1 
and 2 (Table ||), both at z — 3 and z — 0, are in fair 
agreement with these observationally derived quantities. 
For the other experiments, the fraction of stars older than 
z — 3 is significantly lower, which is not only due to a 
low star formation activity at high redshift, but also to a 
higher SFR in the local universe and the effects of cosmic 
variance due to the smaller simulated volumes. 
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Experiment p*(z = 0) p % {z = 3)//?*(0) 



ACDM1 
ACDM2 
ACDM3 
ACDM4 
ACDM5 
SCDM 
BSI 



6.50 
8.60 
9.51 
9.01 
10.9 
12.3 
7.37 



0.252 
0.243 
0.115 
0.116 
0.119 
0.115 
0.041 



Table 3. Present star density (in 10 s h M Q Mpc" 3 ) and 
fraction of stars older than z = 3 for the different numer- 
ical experiments. 



5. Conclusions 

We have performed 66 numerical simulations, grouped in 7 
different experiments, in order to obtain a theoretical pre- 
diction of the cosmic star formation history. We compared 
our predictions with a number of published observational 
estimates. We summarize our results as follows: 



1. 



ACDM shows the best agreement with observations. 
Standard CDM tends to overpredict the cosmic SFR 
density, while BSI fails to form enough stars at high z 
due to its reduced power at small scales. 
Photoionization fixes the onset of star formation in the 
universe, supressing star formation in low-mass objects 
at high z. Supernovae feedbacks are very efficient to 
self-regulate the process of converting gas into stars, 
preventing the most massive halos to form too many 
stars at later epochs. 

Star formation histories are very different for objects 
that end up in the cores of clusters and isolated galax- 
ies. Galaxies in cluster cores have much older stellar 
populations, and their activity is drastically reduced 
when they lose their cold gas reservoirs. 
In contrast to the first observational estimates 



(Madau et al. 1996), but in agreement with most re- 
cent data QSteidel et al. 1999] ), the ACDM model pre- 
dicts almost no drop in the cosmic SFR for 2 < z < 5. 
Star formation seems to be a gradual process with no 
characteristic epoch. 
5. The comoving simulated volumes in our experiments 
are comparable to those covered by observations. Thus, 
one should expect a non negligible error in the observa- 
tional measurements due to cosmic variance associated 
with small volume statistics. This is most important 
when one wants to determine the steep slope of the 
SFR density from z ~ 1 to the present time. 

Selection effects, dust extinction and the conversion 
from luminosities to SFRs are the major sources of uncer- 
tainty in the observational estimates. The small volume 
sampled by the HDF may also introduce some statistical 
bias at high redshift. 

Uncertainty in our numerical experiments comes from 
spatial resolution and the star formation and feedback pre- 
scriptions. The use of Adaptive Mesh Refinement Eulerian 



codes (AMR) will be useful in order to improve res- 
olution in statistically significant cosmological volumes. 
Photoionization and chemical enrichment have been given 
a phenomenological treatment in our code. It would be 
extremely interesting to self-consistently model both pro- 
cesses, computing the intensity of the ionizing background 
at each time-step, as well as metal production and advec- 
tion. 
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